Accurate simulation estimates of cloud points of polydisperse fluids 
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We describe two distinct approaches to obtaining cloud point densities and coexistence properties 
of polydisperse fluid mixtures by Monte Carlo simulation within the grand canonical ensemble. The 
first method determines the chemical potential distribution jit(cr) (with a the polydisperse attribute) 
under the constraint that the ensemble average of the particle density distribution p(a) matches 
a prescribed parent form. Within the region of phase coexistence (delineated by the cloud curve) 
this leads to a distribution of the fluctuating overall particle density n, pin), that necessarily has 
unequal peak weights in order to satisfy a generalized lever rule. A theoretical analysis shows 
that as a consequence, finite-size corrections to estimates of coexistence properties are power laws 
in the system size. The second method assigns j-i(a) such that an equal peak weight criterion is 
satisfied for p(n) for all points within the coexistence region. However, since equal volumes of the 
coexisting phases cannot satisfy the lever rule for the prescribed parent, their relative contributions 
must be weighted appropriately when determining /x(cr). We show how to ascertain the requisite 
weight factor operationally. A theoretical analysis of the second method suggests that it leads to 
finite-size corrections to estimates of coexistence properties which are exponentially small in the 
system size. The scaling predictions for both methods are tested via Monte Carlo simulations of a 
novel polydisperse lattice gas model near its cloud curve, the results showing excellent quantitative 
agreement with the theory. 
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I. INTRODUCTION AND BACKGROUND 

Examples of polydisperse fluids arise throughout soft 
matter science, notably in colloidal dispersions, polymer 
solutions and liquid-crystals. Typically the polydisper- 
sity of such systems is manifest as variation in some at- 
tribute such as particle size, shape or charge, which is 
customarily denoted by a continuous parameter a. The 
state of the system is then quantifiable in terms of a dis- 
tribution p(o~) measuring the number density of particles 
of each a; more precisely, p(a)da is the number density 
of particles in the range a . . . a + da. As such, one can 
regard the system as a mixture of an infinite number of 
particle "species" each labelled by the value of a [1]. 

As has long been appreciated, polydispersity can pro- 
foundly influence the thermodynamical and processing 
properties of complex fluids [2, 3], making a clear eluci- 
dation of its detailed role a matter of both fundamental 
and practical importance. In particular the phase be- 
haviour of polydisperse systems is known to be consid- 
erably richer in both variety and character than that of 
corresponding monodisperse systems (see [4] for a recent 
review). The source of this richness can be traced to 
fractionation effects: at coexistence a polydisperse fluid 
described by some initial "parent" distribution, p(°)(cr), 
may divide into two or more "daughter" phases p( a \a), 
a = 1,2,..., each of which differs in composition from 
the parent. The sole constraint is that the volumetric 
average of the daughter distributions equals the parent 
distribution. 

The occurrence of fractionation can engender dramatic 
alterations to phase diagrams. Insight into the essential 



features of polydisperse phase behaviour can be gained 
by first considering the simpler case of a binary mixture 
of two components whose densities we denote p\ and p2 
(see also Ref. [4]). Let us confine our attention to the case 
of a "dilution line" in the full phase diagram, in which 
we vary (at some fixed temperature) the overall (parent) 
density — pf^ + whilst holding constant the 

ratio of the densities p[ / p\ % i-e. the overall composi- 
tion. These parents thus lie on a straight line through 
the origin in the [p\, /?2)-plane, as shown in fig. 1(a). So- 
called cloud points (marked A and B) delimit the range 
of parent densities for which phase coexistence occurs 
on the dilution line. At cloud point A, the low density 
parent phase coexists with an infinitesimal volume of a 
high density daughter phase A'; while at cloud point B, 
the parent coexists with an infinitesimal volume of a low 
density daughter phase B'. Owing to fractionation, how- 
ever, the compositions p\^ / P2 ( a = 1) 2) of the incipient 
daughter or "shadow" phases differ in general from that 
of the parent: the shadow points (A' and B') lie off the 
dilution line. 

For parents on the dilution line but with densities 
intermediate between the cloud points (e.g. point C of 
fig. 1(b)), two daughter phases (C and C") form. These 
phases occupy finite fractional volumes and their com- 
positions P1/P2 ( a — 1)2) both differ from that of the 
parent. Moreover, the daughter phase compositions vary 
non-trivially as one scans the parent density between 
the cloud points. Consequently, and as a result of the 

lever rule (1— £)pj +£pj 2 ^ = p"p [i — 1, 2), the fractional 
volumes 1 — £ and £ of two phases will in general depend 
non-linearly on the parent density. To fully specify the 
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FIG. 1: Schematic representation of the fractionation be- 
haviour of a binary fluid mixture, as described in the text. 
The parent densities pf , are constrained to the dilution 
line (dashed curve). The curved line indicates the boundary 
of the coexistence region; straight tielines connect coexisting 
phases, (a) At the cloud points A and B, the parent phase 
coexists with the shadow phase A' and B' respectively, (b) 
For a general parent density C in the coexistence region, two 
daughter phases C and C" form. The parent lies on the tieline 
connecting them; this is the geometrical representation of the 
lever rule. 



coexistence properties, one thus needs to determine the 
variation of £ and p^ / (a = 1,2) with n^h 

Turning now to the fully polydisperse case, we con- 
sider a family of parent density distributions (a) = 
n (o) (u) with fixed (normalized) particle size distribu- 
tion f(°\cr) and varying overall number density , the 
value of which parameterizes the location of the system 
on the dilution line. Repeating the above considerations 
for a range of temperatures, one finds that the familiar 
liquid-vapor binodal in the density-temperature plane of 
a monodisperse fluid splits into cloud and shadow curves 
[4], as shown schematically in fig. 2. These mark, respec- 
tively, the density of the onset of phase separation and 
the density of the incipient (shadow) phase [5] . The crit- 
ical point no longer occurs at the maximum of the cloud 
curve but at the intersection of the cloud and shadow 
curves (at which both coexisting phases are identical). 
One interesting implication of this is that even at the 
critical temperature, liquid vapor coexistence can occur 
provided the overall parent density is less than its critical 
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FIG. 2: A schematic fluid-fluid phase diagram for a polydis- 
perse fluid, showing temperature T against density; the cloud 
curve gives the parent density n' - 1 where phase separation 
first occurs while the shadow curve records the density of the 
incipient coexisting phases. 

value. 

In this work we address the problem of accurately de- 
termining phase coexistence properties or, more specif- 
ically, cloud point densities of polydisperse mixtures 
via Monte Carlo (MC) simulation. In a monodisperse 
system, this task is relatively straightforward (see e.g. 
Ref. [6]), because the properties of the coexisting phases 
at a given temperature are the same for all values of the 
overall density lying within the coexistence region (as 
delineated by the binodal). By contrast, for multicom- 
ponent fluids, the occurrence of fractionation implies (as 
we have seen) that the compositions of the phases vary 
across the coexistence region (i.e. as a function of par- 
ent density n^) and one needs to consider carefully the 
repercussions of this for simulation estimates of coexis- 
tence properties. 

The layout of our paper is as follows. In the following 
section we describe the principal computational issues 
arising from fractionation and detail two distinct strate- 
gies for locating cloud point densities within a grand 
canonical ensemble framework. The finite-size scaling 
properties of both methods are analyzed theoretically by 
generalizing to polydisperse systems scaling concepts de- 
veloped originally in the context of monodisperse phase 
equilibria. The predictions for both methods are tested 
in sec. Ill via detailed MC simulations of a novel poly- 
disperse lattice gas model. We conclude in sec. IV by 
comparing and contrasting the relative merits of the two 
approaches. 

II. METHODS FOR CLOUD POINT 
DETERMINATION 

We shall work within the grand canonical ensemble 
(GCE), where the set of chemical potentials /x(cr) is the 
control parameter, while the particle density n and the 
particle size distribution are fluctuating variables. (Tern- 
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perature is assumed held fixed and not written explicitly.) 
The existence of two phases at given chemical potentials 
can be detected from the presence of two separate peaks 
in the probability distribution of the fluctuating order 
parameter, which we take to be the number density n. 
The fractions x^ of probability mass under each 
of the two peaks of p(n) (i.e. the peak weights) are re- 
lated to the pressure difference between the two phases 
as discussed below. For an infinitely large system, where 
the GCE and the canonical ensemble become equivalent, 
a^ 1 ) and x^ would correspond to the fractions of overall 
system volume occupied by each phase, i.e. 1 — £ and £ 
as defined above. 

The density distributions p^ x \a), p( 2 \a) of the two 
phases can be assigned by averaging only over configura- 
tions belonging to either peak of p(n) . These are related 
to the overall average density distribution by 

p(a) =x il) p (1 \a)+x (2) p i2) (a) , (1) 

which can be regarded as a generalized form of the lever 
rule. 

In order to estimate the cloud point for the set of dilu- 
tion line parents p^(a) = f(°\a) , one needs to find 

the value rffi of the parent density at the boundary 
of the coexistence region, where the fractional volume of 
one of the phases, say x^ 2 \ drops to zero. Of course, a 
value of exactly zero is obtained only in the thermody- 
namic limit of infinite system size; the key question is how 
a reliable estimate of can nevertheless be extracted 
from data for finite-sized systems. 



A. Method A 

In method A, proposed originally in [7], we take the 
natural step of tuning the chemical potentials pt(o) such 
that p(a) from (1) equals directly the parent p(°'(a), and 
measure the density-dependence of the peak weight ra- 
tio r = x^/x^\ (We do not detail here the algorithm 
for tuning the p(cr), which is explained in [8]). To un- 
derstand how r will depend on parent density and 
system size, note that r is determined by the difference in 
the grand potential; this is directly related to the pres- 
sure P so that r = cxp(L d (3AP) for large system size 
L. Here d is the dimension of space, (3 = \/k B T and 
AP = _p( 2 ) — PW. The criterion for stable coexistence is 
that r must have a finite value as L — ► oo; the pressure 
difference then has to scale as AP ~ L~ d except in the 
special case r = 1 (see method B below). 

For finite L, metastable coexistence can still be ob- 
served in the density region < where AP = 
0(1), but here r will be exponentially small. To estimate 
n^P from data for a finite system, we use the fact that AP 
is 0(1) and scales linearly with —r$ to leading order 
near the cloud point, and hence lnr <~ L d (n^ — n^P). 



This applies for < , while above n£j one has 
lnr = 0(1). Thus the derivative (d/dn^)\nr should 
drop from an 0(L d ) plateau to O(l) around 
In the second derivative — (d/dn^) 2 lnr this drop will 
manifest itself as a peak, whose position serves as an es- 
timate for n^p. We next derive the finite-size scaling of 
the location and shape of this peak. 

Our starting point is the rigorous result of [9] for first- 
order phase transitions driven by some field h. The au- 
thors of [9] showed that, if the transition is at h = 0, the 
free energy density around this point is given by 

f(h) = -L- d k B T\n (e- Ld M W W + e^WW) (2 ) 

up to exponentially small corrections of order 
L~ d exp(— const x L), which can be neglected. The 
function f^\h) (f^(h)) is the thermodynamic free 
energy of the first (second) phase in the regime where 
that phase is stable; elsewhere, i.e. for values of h 
where the phase would be only metastable in an infinite 
system, it can be chosen as the continuation with smooth 
derivatives (up to at least 3rd order [9]) of the stable 
free energy. Intuitively, the approximation (2) tells us 
that the partition function close to the phase transition 
can be obtained by adding the partition functions of the 
coexisting (stable or metastable) phases. Formally, it is 
valid only for values of \h\ <C L^ 1 ; but in fact we will be 
interested in rather smaller values of h <~ L~ d InL where 
the smaller phase is not yet exponentially suppressed, so 
that this is not a restriction. 

For a polydisperse system within the GCE, the ana- 
logue of the field h is the set of chemical potentials. Con- 
ceptually it is easiest to think first of particle sizes dis- 
cretized into a large but finite number M of bins; there 
are then M chemical potentials and densities which we 
write as vectors p and p, respectively. More precisely, 
we take p. as the difference of the chemical potentials 
from those at the desired cloud point, so that the latter 
is located at p = 0. Assuming that the result (2) can be 
extended to situations with M field variables instead of a 
single one we have then for the negative grand potential 
density, which is nothing but the pressure, 

P{p) = L- d k B Tln (e Ld f* pW M + e L ^ pW ^) (3) 

The single-phase pressures expanded to second order in 
p, read (a = 1,2) 

p^(p) = p cl + p i ; l ) -M + ^-x (a V (4) 

Here P c i is the coexistence pressure at the cloud point, 
which is common to both phases, while p^ and p^ 
are the (vectors of) particle densities at the cloud point; 
p^ is then equal to the parent density vector p^? at 

(2) 

the cloud point, while p y cl ' is the shadow density vector. 
The matrices x (a) = V M p (a) = V M V M P (a) are the sus- 
ceptibilities of the particle densities to chemical potential 
changes, again evaluated at the cloud point. 
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Taking the chemical potential derivative of (3) we have 
for the overall density vector 



P = 



(5) 

This is of the form (1) with the natural correspondences 
p (a) (p) = V„P(°>(/i), = 1/(1 + r), x& = r/(l + r) 
and lnr = L d (3[P( 2 \p) — P^ l \p)}. We now expand near 
the transition, keeping terms to the same order as in (4). 

With the abbreviations Ap c i = p^' — p^' and Ax = 
X^~X^ an( i using that method A imposes p = n^f^, 
where f (°) is the normalized parent size distribution, one 
has 



l (0)f(0) 



■x ( V 



1 + r 

lnr = L d f3(Ap ■ p + p ■ A X p) 



(Ap cl + A X p) (6) 
(7) 



After eliminating p, these relations determine the de- 
pendence of r on that we seek. To make progress, 
we anticipate that the peak in the second derivative 
— (d/dn^) 2 lnr will occur at a point where r = 0(L~ d ). 
From (7) this implies that p = 0(L~ d \nL) is also small. 
Keeping only leading order terms in the small quanti- 
ties r and p and using that p^ 

) =n (0>f(0) 

, our previous 

relations then become 

(„(0) _ „(°))f(0) 



X {1) P + rAp ci 



lnr = L d f3Ap cX ■ p 



(8) 
(9) 



Solving the first for p and inserting into the second then 
gives 

lnr = L d (3Ap cl • ( X (1) ) _1 [(n^ - n^)f (0) - rAp cl ] (10) 

To absorb the numerical coefficients and make the parent 
density dimensionless we define 



n<°> 



arL d 

bL d (nW -n^) + ln(aL d ) 



with 



a = /3Ap cl • (x (1) ) _1 Apci 
b = /3A Pcl • (x (1) )" 1 f (0) 

The relation (10) then becomes just 
n (0) = z + In z 



(11) 
(12) 



(13) 
(14) 



(15) 



Differentiating w.r.t. n(°) gives (d/dh^) Inz = (z + 1) 1 
and -(d/dn^) 2 \nz = (z + \y 2 z (8/ On™) In z = z(z + 
l) -3 . Bearing in mind that lnz and lnr differ only by a 
constant, we therefore arrive at a universal large-L scal- 
ing form for our second-derivative plot 
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lnr 



(1 + z) 3 ' 



z + lnz (16) 



which is parameterized by z. All dependence on system 
details is encoded in the two numerical constants in the 
definition (12) of the scaled parent density. The curve 
(16) has its peak at z = 1/2 so that, from (11), r is 
0{L~ d ) in the region of interest as anticipated in our 
derivation. The position of the peak on the horizontal 
axis is n(°) = (1/2) - In 2. The scaling (12) then im- 
plies that the cloud point estimated from the peak posi- 
tion has finite-size corrections of order L~ d In L, while the 
peak width and height scale as L~ d and L 2d , respectively. 
We will find these scalings, and indeed the full shape of 
the master curve (16), confirmed in the simulation data 
shown below. 



B. Method B 

While the scaling analysis described above provides a 
detailed picture of the finite-size corrections that arise 
when using method A to estimate the location of the 
cloud point, it would clearly be desirable from a prac- 
tical point of view to reduce these corrections and ide- 
ally make them exponentially small in system size. For 
monodisperse phase coexistence, this is achieved by mea- 
suring the densities of the coexisting phases at the spe- 
cial point where the peaks of p(n) have equal weights 
X W = x ( 2 ) = 1/2, i.e. r = 1. At this point the pressure 
difference vanishes. In method A, on the other hand, 
the pressure difference is AP = L~ d kBTlnr <~ L~ d \nL, 
and this causes the relatively large finite-size corrections. 
This observation suggests that one should also consider 
coexisting phases with r = 1 in the polydisperse case. Of 
course, one can then no longer require the parent density 
distribution to equal the overall density distribution in 
the system, since the latter is an equal mixture of (a) 
and p( 2 \a). One has to allow more generally 



p<°>(a) = (l-Op (1 V)+fr (2 V) 



(17) 



where £ is a parameter to be determined; the cloud point 
is estimated as the parent density at which £ reaches zero. 
The results do, however, have physical meaning also for 
other values of £ in the range < £ < 1: they then 
provide estimates of properties of the coexisting phases 



for parent densities 



,(o) 



/ithin the coexistence region, 



with £ estimating the fractional volume of phase (2). In 

simple cases, the cloud point n^P could in fact be esti- 
mated by linearly extrapolating a few measurements of 
£(n(°)) to £ = 0. In a monodisperse system this would be 
exact since £(n.(°)) varies strictly linearly across the co- 
existence region. In the polydisperse case, on the other 
hand, C(n (0) ) is a nonlinear function and can exhibit very 
significant curvature near the cloud point [4, 7, 10]. Lin- 
ear extrapolation is then unreliable and best avoided in 
favour of direct determination of the density where 
£ = 0, as explained above. 

We will call this approach "method B". In practice, 
it is implemented as follows. A parent density is fixed, 
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along with a trial value of £. The chemical potentials are 
then tuned until the density distributions of the two co- 
existing phases satisfy (17). One measures r = x^ /x^; 
if this deviates from r = 1, £ is adapted (e.g. using a bi- 
section method) and the process is iterated until r = 1 to 
numerical accuracy. If a solution with < £ < 1 is found, 
we are in the coexistence region. The parent density is 
then reduced and the process repeated until £ drops to 
zero or no solution with positive £ can be found. 

It is straightforward to adapt the above scaling analysis 
to show that, within method B, finite-size corrections 
are indeed exponentially small. The condition (17) with 
£ = together with r = 1 gives as the analogues of (8,9) 



(n 



(o) 



)f(0) 



= L d (3Ap-n 



(18) 
(19) 



Solving the first and inserting into the second gives 

0=(n(°)-n(° ) )L d /3A Pcl (x (1) )- 1 f (0) (20) 

which (barring accidental vanishing of the constant fac- 
tor, equal to b from (14)) is satisfied only for 
Finite-size corrections therefore arise only from terms 
that we had discarded from the outset in our analysis; 
these are exponentially small. Note that in a practical 
implementation it is important that £ is determined to 
high accuracy; indeed, the analysis above assumes that 
there is no error in the value of £. It is easy to see that, 
if instead £ was found only with an accuracy of 0(L~ d ), 
then finite-size corrections of the same order as in method 
A arise. 

To summarize, method B is algorithmically a little 
more involved than method A because it requires for each 
parent density an "inner loop" over £, but compensates 
for this by producing much smaller finite-size corrections. 
This it achieves by forcing coexisting phases to have iden- 
tical pressures. The parameter £ has to be introduced to 
ensure that the parent distribution is still obtained as an 
uneven mixture of the two phases. Note that this is a pe- 
culiarity of the polydisperse case: no such parameter is 
necessary for monodisperse systems since the properties 
of the coexisting phases are the same everywhere within 
the coexistence region. In a polydisperse scenario, on the 
other hand, the coexisting phases do change [4], and so 
it is important that the mixing proportions appropriate 
to the chosen parent density are maintained. 



III. APPLICATION TO A POLYDISPERSE 
LATTICE GAS MODEL 

A. Model definition 



choice of a lattice-based rather than a continuum model 
was made on the ground of computational tractability: it 
permits the study of a larger range of system sizes than 
is feasible for continuum models. We do not expect the 
scaling predictions for the cloud point estimates to be 
affected by the presence, or otherwise, of a lattice. 

Our polydisperse lattice gas (PLG) model is defined 
within the grand canonical ensemble by the hamiltonian: 

H = - <ro / c i (iT)c J (o')-Y i vi(tT)c i {<T). (21) 

Here a is the particle "species" , whose chemical potential 
is jU(<r), while Ci(a) is the number of particles of species 
a at site i, for which we impose a hard-core constraint 
such that J2a c i( <J ) = or 1. The instantaneous density 
distribution follows as p(a) = L^ d J2i c i( a )j with d = 3 
in the simulations reported below; i runs over the sites 
of a periodic lattice i = 1,...,L , assumed cubic in this 
work. The sum in the first term on the right hand side of 
(21) similarly runs over all pairs i,j of nearest neighbor 
sites, as well as over all combinations of a and a'. 

For a study of this model under conditions of fixed 
polydispersity, one requires knowledge of the chemical 
potential distribution /j,(<j) corresponding to some pre- 
scribed form of the ensemble averaged density distribu- 
tion p{a). In the context of method A, one tunes p(<r) 
such that p(a) = (a) , while for method B one simul- 
taneously tunes p(cr) together with the parameter £ to 
satisfy (17) and the equal peak weight criterion for p{n). 
Such tuning can be efficiently achieved by a combina- 
tion of a non-equilibrium Monte Carlo procedure and his- 
togram extrapolation techniques, as has been described 
previously elsewhere [8] . 

We have studied the dilution line properties for a par- 
ent distribution having the Schulz form: 



/ (0 V) 



1 



z + i 



z+i 



o z exp 



Z + l 



(22) 

The mean value of the distribution a defines the unit 
length, while the parameter Z controls the width of the 
distribution. We fixed the latter to be Z — 50, resulting 
in a dimcnsionless degree of polydispersity 



5= 



af 



1 



^/ZTT 



~ 14% 



(23) 



Lower and upper cutoffs were imposed on the distribution 
at a = 0.5 and 1.4 respectively, and the distribution was 
normalized accordingly. 



In order to test the predictions for the finite-size scaling 
properties of the two methods detailed above, we have 
performed a systematic Monte Carlo simulation study of 
a novel lattice gas model for a polydisperse fluid. The 



B. Simulation results 

A determination of the critical point parameters for 
the PLG model using well-established finite-size seal- 
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ing methods [6] found the critical temperature to be lo- 
cated at (n? } ,T c ) = (0.521(1), 1.171(1)) in reduced units. 
This is to be compared with the critical parameters of 
the monodisperse (Ising) lattice gas (0.5,1.127955) [11]. 
Thus the inclusion of polydispcrsity of the form (21) is 
seen to raise both the critical temperature and the criti- 
cal density. Moreover it splits the liquid-gas binodal into 
well separated cloud and shadow curves (cf. fig. 2) in 
a manner similar to that occurring in continuum fluid 
models for which polydispersity affects the interparticle 
interaction strength as well as its range [7, 12]. As a con- 
sequence, the critical point lies below the maximum of 
the cloud curve and phase coexistence can be observed 
even at T = T c , provided that < 

In view of this we have adopted the critical temper- 
ature as a convenient reference point and have studied 
coexistence along the dilution line p^(cr) = n^f^°\a) 
for fixed T = T c . This involved performing a series of 
MC runs starting from the critical density and reducing 
n(°) in a stepwise fashion. Multiple histogram reweight- 
ing techniques were employed to estimate the chemical 
potential distribution at each successive step, while use 
of multicanonical preweighting techniques [13] ensured 
that both coexisting phases were efficiently sampled in 
the course of each simulation run (see ref. [12] for a fuller 
account of this procedure) . 

The dilution line was scanned in this manner for lat- 
tices of sizes L = 10, 12, 15, 18, 21 using both methods 
A and B. The GCE simulations directly yield the form 
of p(n) corresponding to either case. For method A, one 
observes (figs. 3(a) and 3(a')) that as is reduced from 
its critical point value, the peaks in p(n) separate, while 
the valley between them deepens. This is accompanied 
by a gradual transfer of weight from the liquid to the gas 
peak. We extract the ratio of the peak areas at a given 
n (o) f rom p(n) via 
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Here n* is a convenient threshold density intermediate 
between vapor and liquid densities, which we take to be 
the location of the minimum in p(n) . 

The form of p(n) obtained using method B is shown in 
fig. 3(b) for a selection of values of n^°\ Here, by virtue 
of the appropriate tuning of the parameter £, the vapor 
and the liquid peaks maintain equal weights throughout 
the coexistence region. The estimate of the ratio of frac- 
tional phase volumes is obtained simply as = £/(l — £). 
We use the subscript £ here to distinguish this quantity 
from the probability mass ratio r as extracted from (24), 
the latter being more directly related to the pressure dif- 
ference between the phases as explained above. 

In fig. 4, we plot the dependence of r on parent density 
n(°) as obtained from method A for the various system 
sizes we have studied, alongside the results for from 
method B. The curves for method A (thick lines) show a 



FIG. 3: (a) Estimates of the form of p(n) for a selection 
of values of n' ' within the liquid-vapor coexistence region 
as obtained from the GCE simulation data obtained accord- 
ing to method A. The associated estimates of r are, in order 
of decreasing n (0) : r = 0.4483, 0.2682, 0.1643, 0.0901, 0.0379. 
(a') The distributions for a selection of small values of n^ ' 
displayed on a logarithmic scale, (b) The corresponding esti- 
mates of p(n) obtained according to method B (for which 
the peak weight ratio is constrained to unity). The as- 
sociated estimates of are, in order of decreasing n' - 1 : 
r ( = 0.4158,0.2411,0.1477,0.0822,0.0335. All the plots re- 
fer to system size L = 15. 



strong i-dcpcndcncc within the metastable coexistence 
region which borders the cloud point at small . As 
n(°) is increased, however, and the cloud point is ap- 
proached, the curves cross over to their L-indcpendent 
limit values in the coexistence region. 

Looking in more detail at the finite-size scaling proper- 
ties of the curves for r(n^ ) from method A, the analysis 
of sec. II A shows that the cloud point region is associated 
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FIG. 4: The value of r as estimated from method A (la- 
belled thick curves) and of from method B (unlabelled, 
thin curves); both are plotted against parent density nf-°\ 
The inset shows in greater detail the results for method B in 
the vicinity of the cloud point density. An estimate of the 
statistical errors of the procedure were obtained via a block 
analysis. The errors for the smallest and largest system size 
are indicative of the statistical uncertainty of our results. 
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FIG. 5: Points show the estimated second derivative 
— (d/dn^) 2 lnr obtained within method A. The solid lines 
superimposed on these data points are scaled fits of the form 
(16). Insets show the //-dependence of the peak width and 
height plotted in terms of the predicted scaling variable and 
accompanied by a least squares fit. 



with a peak in the second derivative ~{d/dn^) 2 lnr. In 
Fig. 5 we plot this quantity for the various systems sizes 
studied. Superimposed on each plot is a suitably scaled 
form of the predicted universal master curve, (16). The 
non-universal scale factors for the height and the width of 
the master curves implicit in the fitting shown are plot- 
ted in the insets against the predicted scaling variables 
L 2d and L~ d respectively (cf. sec. II A). Clearly there 
is good agreement with the predictions for both the gen- 
eral shape of the peak and the scaling of its width and 



height. Some discrepancies between the observed and 
measured peak shape are apparent on the low density 
side, well away from the peak maximum, particularly 
for the smaller system sizes. These are presumably at- 
tributable to the breakdown of the validity of the linear 
expansion in the density difference used in the derivation 
of (16). 

The estimates of r^(n^) deriving from application of 
method B, as shown in Fig. 4, exhibit a qualitatively 
different behavior. As the density is increased from the 
cloud point, one observes a rapid, near- vertical increase 
in In rj from a large negative value (where is essentially 
zero) to values of 0(1). There is only a weak dependence 
of this behaviour on the system size (see inset of fig. 4). 
Within this method the cloud point density can thus be 
directly read off as the lowest density at which an equal- 
peak-weight solution exists to the numerical procedure 
described in Sec. II B. We note that, as expected, r from 
method A and from method B converge to similar 
values for parent densities well within the coexistence 
region. 

Fig. 6 compares for the two methods the finite-size be- 
haviour of the cloud point density estimates. In the case 
of method A, the estimate for a given L derives from the 
position of the peak in the second derivative plot (fig. 5). 
Fig. 6 confirms that (as predicted) these estimates devi- 
ate from their limiting value by a correction 0(L~ d In L). 

A least-square fit yields nf^ = 0.20266(12) as the best 
estimate of the cloud point density (see also inset). In the 
case of method B, the cloud point is estimated as the low- 
est value of for which an equal-peak-weight solution 
for p{ri) can be found. The finite-size corrections to this 
estimate are expected to be exponentially small in the 
system size and our results (fig. 6(b)) are consistent with 
this, yielding a cloud point estimate = 0.2028(1). 
Indeed the corrections appear to be so small that even 
for the smallest system size studied (L = 10) the esti- 
mate of the cloud point density obtained using method 
B deviates from the limiting value by just 0.5%. This 
compares with a relative deviation of approximately 10% 
for method A. 



IV. CONCLUSIONS 

We summarize. We have presented and tested two 
distinct finite-size scaling strategies for estimating coex- 
istence properties and cloud points in polydisperse flu- 
ids within a grand canonical ensemble framework. Both 
are, of course, also directly applicable to multicompo- 
nent mixtures of many discrete species. The first ap- 
proach, "method A" takes the natural step of constrain- 
ing the ensemble average of the density distribution p(a) 
to match the prescribed parent form. However, in order 
to satisfy the lever rule, this necessarily leads to unequal 
peak weights in the order parameter distribution function 
p{n). For systems of finite size, this translates to a finite 
pressure difference between the coexisting phases, which 
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in turn engenders finite-size corrections to estimates of 
the cloud point density which are powers of the system 
size. 

The second approach, "method B" strictly imposes 
equal peak weights for p(n) , and hence pressure equality 
in systems of finite size. The coexistence properties of the 
desired parent are obtained by appropriately weighting 
the relative contributions of the coexisting phases to the 
overall density distribution in such a way that the lever 
rule is satisfied. We have detailed how to determine this 
weight factor operationally and shown that method B has 
finite-size corrections to the coexistence properties which 
are exponentially small in the system size. As such, and 
notwithstanding the need to determine the phase weight- 
ing factor £, method B is clearly superior to method A. 
It can be regarded as a generalization to multicomponent 
mixtures of the well known equal peak weight criterion 
[9] developed in the context of monodisperse systems by 
Borgs and coworkers. 
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FIG. 6: Finite-size scaling of the cloud point estimates deriv- 
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gion of large L. (b) The estimates deriving from method B, 
plotted against L~ d exp(-nL) with k = 8.988 x 10~ 3 . Ex- 
cept where error bars are shown, statistical uncertainties are 
smaller than the symbol sizes. 
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